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ABSTRACT 



We compare the results for the dynamical evolution of star clusters derived from 
anisotropic gaseous models with the data from A^-body simulations of isolated and one- 
component systems, each having modest number of stars. The statistical quality of A^- 
body data was improved by averaging results from many A^-body runs, each with the 
same initial parameters but with different sequences of random numbers used to initialize 
positions and velocities of the particles. We study the development of anisotropy, the 
spatial evolution and energy generation by three-body binaries and its A^-dependence. We 
estimate the following free parameters of anisotropic gaseous models: the time scale for 
coUisional anisotropy decay and the coefficient in the formulae for energy generation by 
three-body binaries. To achieve a fair agreement between N-body and gaseous models for 
the core in pre- as well as in post-collapse only the energy generation by binaries had to 
be varied by N. We find that anisotropy has considerable influence on the spatial structure 
of the cluster particularly for the intermediate and outer regions. 
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1 INTRODUCTION 



One of the grand challenges of theoretical astrophysics is to understand the dynamics 
of globular star clusters. Aside from the intrinsic interest of understanding the behaviour of 
large A'"-body systems, the importance of the problem stems from its relation with current 
research into the stellar content of globular star clusters. (See, for example, many of the 
papers in Janes 1991). Unfortunately the direct simulation of such rich stellar systems 
with A'"-body modelling is not yet feasible (Hut, Makino & McMillan 1988). The gap 
between the largest useful computer models {N < 10000) and the median globular star 
cluster (N ~ 5 x 10^) can only be bridged at present by use of theory. There are two 
main classes of theory (cf. Saslaw 1985): (i) Fokker-Planck models, which are based on 
the Boltzmann equation of the kinetic theory of gases, and (ii) gas models, which can be 
thought of as a set of moment equations of the Fokker-Planck model. 

These simplified models are the only detailed models which are directly applicable 
to large systems such as globular star clusters. But their simplicity stems from many 
approximations and assumptions which are required in their formulation, and it is not 
clear how well these correspond to the real systems in nature. The usual assumption of 
spherical symmetry contradicts the asymmetry of galactic tidal fields. Almost all Fokker- 
Planck calculations assume isotropy of the velocity distribution, which contradicts the 
evidence of proper motions and of model-building. The treatment of stellar escape in 
Fokker-Planck models has long been problematic, and it is never considered in gas models. 
Three- and four-body encounters can only be modeled within our limited knowledge of the 
relevant scattering cross-sections. 

In view of all these uncertainties, the reliance which can be placed on these simplified 
models is doubtful. For this reason it is of importance to test the predictions of these 
models against results which are not subject to this wide class of simplifying assumptions 
and approximations, i.e. A^-body models. Despite the central role of the simplified models 
in much recent research on cluster dynamics, little has been done to test their validity. 
There is the already classical comparison between the fluid-dynamical model of Larson 
(1970), Monte-Carlo simulations, and direct A-body integrations of 100 and 250 particles 
(Aarseth, Henon & Wielen 1974, see also Aarseth & Lecar 1975), which were at that time 
maximal particle numbers from the viewpoint of computational resources. 

In the following years there has been much further improvement of the theoretical mod- 
els as well as of computational resources. Direct numerical solutions of the Fokker-Planck 
equation substituted the Montc-Carlo models, and gaseous models with heat flux closure 
improved Larson's original fiuid-dynamical approach (Lynden-Bell & Eggleton 1980, Heg- 
gie 1984). Bettwieser & Sugimoto (1985) tried to check the main assumption of the gaseous 
model, the heat conductivity, by using a direct N — 1000 model. Their A-body results, 
although in fair agreement with expectations in pre-coUapse, suffered from large statistical 
fluctuations especially in post-collapse due to the still too small particle number. 

With the advent of several large parallel computing facilities as well as fast vector 
computers in general it is now possible to extend direct A-body calculations in two respects: 
first to improve the statistics for rather low A (up to a few thousand) by computing many 
independent models simultaneously, one on each processor, improving the statistics then by 
averaging. On the other hand the use of one of the fastest available vector supercomputers 
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(CRAY YMP) yielded for the first time high accuracy A'"-body models for as much as 
10000 particles throughout most of the core collapse phase (Spurzem & Aarseth 1993). 

This paper is designed to complement a larger survey of star cluster evolution (Giersz 
& Heggie 1993ab), which compares A^-body models with isotropic gaseous and Fokker- 
Planck results, by one particular aspect; this is to check the reliability of a generalized 
gaseous model which includes the possible generation of anisotropy in the radial and tan- 
gential dispersions of the star cluster. On the other side there is a feedback from the 
A-body models back to some adjustments in the gaseous model code in order to reach 
optimal agreement of both models. 

In the following section theory and numerical solutions of anisotropic gaseous models 
are elaborated. Sect. 3 presents some information on the A-body models used, to an extent 
that is necessary here for understanding this paper; more details will be given elsewhere 
(Giersz & Heggie 1993a). Sect. 4 gives an account of the results, which are finally discussed 
and complemented by concluding remarks in the final section. 



2 ANISOTROPIC GASEOUS MODELS OF STAR CLUSTERS 
2.1 The model 

Observational fits of globular clusters (cf. e.g. Lupton & Gunn 1987, Lupton, Gunn 
Sz Griffin 1987) and direct A-body calculations show that there is a considerable amount 
of anisotropy (cr^ > cr^) in their halo. Including the effects of anisotropy into models 
of the dynamical cluster evolution has posed some difficulties in the past. The only 2-D 
orbit average Fokker-Planck models for anisotropic single mass clusters with and without 
a massive central object (Cohn & Kulsrud 1978, Cohn 1979, Cohn 1985) have not been 
used further and there were a variety of anisotropic gaseous models with different closures 
(Bettwieser 1983, Bettwieser & Spurzem 1986); however their relation to real many-body 
systems and the quality of the numerical solutions remained unclear. Recently, however, 
self-similar anisotropic models for regular pre-coUapse as well as for singular post-collapse 
clusters were obtained (Louis & Spurzem 1991, henceforth LS). This is an occasion to 
revisit the quality of the numerical solutions of the anisotropic gaseous model and to 
utilize such a model further to compare its predictions with the results of direct A-body 
calculations, with emphasis on the results concerning the anisotropy. 
For the sake of completeness and to clarify small differences in the model and notation 
compared with LS we will present here the full set of equations used; dependent variables 
are the mass Mj. contained in a sphere of radius r, the local mass density p, radial and 
tangential pressure PriPti bulk mass transport velocity tt, and transport velocities Vr, Vf of 
the radial and tangential energy, respectively. As auxiliary quantities we use the radial and 
tangential 1-D velocity dispersions = Pr/ cr^ = pt/p, the average velocity dispersion 
(T^ = (cT^ -I- 2(jf )/3, the anisotropy A — 2 — 2af/a'^ (note that A — 6a/(l -|- 2a), where a is 
the anisotropy measure used in LS), and the relaxation time 
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in the definition of Larson (1970), where N is the total particle number of the star cluster, 
m the individual stellar mass and 7 a numerical constant whose value will be discussed 
below. The equations are 

^=4..V (2) 
%^^l(>-r^) = (3) 



?^+^?}i + ^ + l?EL + 2f^ = (4) 
ot or p or pr 



dt dr dr r'^ dr r 3 Aa^a \ / bins 

(5) 

Vr -U-\ = (7) 

47TGpT dr ^ ' 

Vr = Vt (8) 

They are nearly equivalent to model A of LS (equal velocities for the transport of radial 
and tangential thermal energy) . Note that model A (the so-called "1 flux" model) is very 
similar to the gaseous model of Bettwieser & Spurzem (1986), whereas model B of LS ("2 
flux" model) is related to that of Bettwieser (1983). The net transport velocities for radial 
and tangential energy {vr — u) and {vt — u) can be derived from the energy fluxes Fj. and 
Ff used as variables in these papers by dividing out a convenient multiple of the relevant 
pressure {2pt for {vt — u), 3pr for {vj. — u) ). The reader interested in more details about 
this and the connection of the variables to moments of the stellar velocity distribution is 
referred to LS. 

Apart from the energy generation term due to hard binaries discussed below there is one 
more difference to LS in keeping the hydro dynamical terms in Eq. (4); for all applications 
comparing with A/"-body calculations this did not make any differences, however the sta- 
bility of the code at very small timesteps and high central densities (near core bounce) is 
enhanced; if one follows core collapse over many orders of magnitude increase in central 
density these terms slightly affect the value of the anisotropy in comparison with the self- 
similar models of LS. Hence they were omitted only for test calculations to be compared 
quantitatively with the results of that paper. 

The numerical constants Aa and A occurring in Eqs. (5) to (7) are related to the 
timescales of collisional anisotropy decay and heat transport, respectively. A is related to 
the standard C constant in isotropic gaseous models (see e.g. Heggie & Stephenson 1988) 

by 

A = (8) 
10 ^ ^ 
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Ta is the anisotropy decay timescale for an anisotropic local velocity distribution function; 
in the Appendix it is outlined how one gets the final result Ta = lOT/9, provided a 
particular velocity distribution function is assumed. Since the real distribution function is 
not completely known within the framework of a gaseous model a free numerical constant 
Aa is introduced in Eqs. (5) and (6); the values chosen for Aa and A will be discussed 
in comparison with the direct A^-body calculations later. Note that in LS and all other 
previous papers the timescale Ta for an isotropic stellar background distribution as derived 
by Larson (1970) was used, which is 5T/6. Thus all models of LS correspond to a choice 
of Aa = 3/4 here. 

For the following comparison A/"-body models consisting of equal gravitating point masses 
are chosen, where any effects of finite size stars or stellar evolution are neglected. In that 
case the dominant energy source, which finally halts core collapse is the energy generated 
by formation and hardening of three-body binaries, which we take in accord with the 
standard ansatz (Goodman 1987) as 



This is an isotropic energy input. A 6'-function has been attached here to illustrate the 
possibility to start the binary energy generation not from the very beginning but at a time 
tbo, which is close but not identical to the time of core collapse and will be elaborated in 



Although the standard ansatz chooses Cb = 90, we keep this value here as a free parameter 
for reasons which will also be discussed in Sect. 4. Such a type of energy source to model the 
average energy generation by three-body binary generation and hardening was first used 
by Bettwieser & Sugimoto (1984) and Heggie (1984) for their isotropic gaseous models. 
Note that our anisotropic gaseous model equations are identical up to second order to 
moment equations of the Boltzmann equation with a Fokker-Planck collisional term. The 
closure equations Eqs. (7) and (8) are related to the analogous equation of heat transfer 
in a gas, except that the timescale occurring in the conductivity here is the appropriate 
stellar dynamical relaxation time (Lynden-Bell & Eggleton 1980); thus such models are 
denoted as gaseous models to distinguish them from higher order fiuid dynamical models 
(e.g. Louis 1990). 

Altogether the anisotropic gaseous model itself contains three parameters. A, 7, Aa; in- 
cluding the binary energy generation adds another two parameters, namely Cf, and tbo- 
It is shown in a separate paper comparing isotropic gaseous models with direct A'^-body 
calculations (Giersz & Heggie 1993a) that for all models (including also direct solutions of 
the Fokker-Planck equation) values close to 7 = 0.11 and C = 0.104 (i.e. A = 0.4977), give 
the best agreement; the isotropic and anisotropic gaseous models differ only marginally 
in the evolution of the inner mass shells during pre- and post-collapse evolution; this is 
shown in Fig. 1, where an isotropic model obtained with Heggie's code is compared with 
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an anisotropic 1 flux model of Spurzem's code. Since the behaviour of the inner Lagrangian 
radii determines the best values for 7 and C as discussed by Giersz & Heggie (1993a) we 
use here just the same values. 

It is the focus of this paper to show how a comparison with direct A^-body calculations 
enables us to find unambiguously an optimal set of the remaining parameters Aa, C*;,, and 
tjjQ in order to achieve fair agreement with the A^-body results. 

Without binary energy generation the quasi-static evolution of a system governed by Eqs. 
(2) to (8) depends only on the product A- Aa; different choices for the total particle number 
result merely in a rescaling of time. With binary energy generation, however, increasing 
means a smaller binary energy generation; thus the post-collapse evolution is qualitatively 
different for varying particle number, ranging from simply steady reexpansion for low 
N , through periodic oscillations at intermediate N and probably chaotic large-amplitude 
gravothermal oscillations (originally detected by the isotropic gaseous model of Bettwieser 
& Sugimoto, 1984) at very large N (Heggie & Ramamani 1989, Cohn, Hut & Wise 1989). 

2.2 Numerical Solution of the Anisotropic Gaseous Model Equations 

For the numerical solution of the model equations and comparisons with direct N- 
body results standard A'"-body units were used, where G = 1, the total mass of the system 
M = 1, and the total energy of Plummer's model, which was used as initial model, is 
E = —1/4; in these units the scaling radius of Plummer's model is a = 37r/16. The 
equations were discretized on an Eulerian mesh with 200 logarithmically equidistant grid 
points (for some test runs 400); they were distributed between a minimal and maximal 
radius of 2.06 • 10~^ and 144 in the above units. 

The equations used for discretization and numerical solution are equivalent but not iden- 
tical to those of Eqs. (2) to (8). First instead of the physical quantities we used logM^., 
logp, \ogpr and iogpt as well as the net energy transport velocities Vr — u and vt — u 
together with u as dependent variables; in the Eulerian mesh scheme r is the indepen- 
dent variable. In order to achieve the utmost accuracy of mass and energy conservation 
a partially implicit scheme was used where the equations were formulated in terms of 
(^X + (1 — C)^*7 where X is the actual value and X* the one of the previous timestep; 
such a scheme is numerically unstable with the optimal value of C = 0.5, thus ( — 0.55 was 
used. The difference equations were discretized on a logarithmically equidistant mesh and 
solved iteratively with a Newton-Raphson-Henyey method. For convenient discretization 
in the logarithmic variables the divergence terms occurring in Eqs. (3), (5), and (6) were 
all split up. As a criterion for choosing the timestep we used that the positive logarithmic 
quantities should not change by more than 0.05, complemented by the condition that the 
timestep should not exceed the central relaxation time, which ensures good accuracy at 
the turning points from pre- to post-coUapsc. As boundary conditions we imposed time- 
independent logarithmic gradients of density and pressures at the outer, linear variation of 
the velocities with respect to r and of the mass with respect to at the inner boundary. 
The conservation of energy and mass in all results presented in the following section, 
related to the total values of 0.25 and 1, respectively, was always smaller then or equal to 
0.25%; although the strict mass conservation of earlier codes was lost (due to splitting the 
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divergence term in Eq. 3 and the different boundary condition) tlie mass error remained 
small as stated; therewith, however, a much better quality of energy conservation was 
achieved than in previously published anisotropic models of Bettwieser & Spurzcm (1986) 
and Bettwieser (1983). Although they discretized the divergence terms in Eqs. (5) and 
(6) in a strictly conservative manner this did not ensure that the total energy, consisting 
of thermal and gravitational energy (and bulk kinetic energy, which is here not important) 
is strictly conserved; the energy conservation becomes worse due to the non-conservative 
form of terms if one has separate energy equations for radial (Eq. 5) and tangential energy 
(Eq. 6). On the contrary the logarithmic variables improve energy conservation since they 
vary linearly whenever there is a power law, which turned out to be the more important 
effect. 

In fact it occurred that the achieved high level of accuracy was necessary here, because for 
the previous models the errors were larger than the remaining deviations between gaseous 
model and direct A^-body calculation. 

The numerical method was tested by comparing with the known self-similar solutions of 
LS; for that purpose the spatial resolution was enhanced by choosing 400 grid points with 
an innermost radius of 1.11 ■ 10~^ (to be comparable with LS the two first hydrodynamic 
terms in Eq. (4) had been omitted for these calculations only). A pure pre-coUapse solution 
was followed (without binary energy generation) over an increase in central density of 13 
orders of magnitude; the energy error here remained below 1% until the density increased 
roughly 8 orders of magnitude - thereafter the resolution of the core became poorer and 
consequently the energy error increased further. It would be possible to increase the 
accuracy even after such a large growth of central density by either more meshpoints or 
moving them with e.g. the shrinking core radius. This is not useful for our purpose here, 
since we do not follow the system deep into core collapse. 

Figs. 2ab show the radial profiles of the logarithmic density power law index a and of 
the anisotropy A for a number of models with different evolutionary time; the onset of 
the self-similar collapse phase is clearly visible. In Figs. 3ab the nearly constant values of 
a and A after the self-similar phase has been entered are shown as a function of AaA in 
comparison to their values in the self-similar anisotropic pre-coUapse models of LS, to the 
isotropic gaseous models of Lynden-Bell & Eggleton (1980), and to the direct solution of 
the isotropic Fokker-Planck equation (Cohn 1980). (Since LS only published two different 
A parameters, additional results for different AaA have been obtained by P.D. Louis, priv. 
communication). Whereas in LS A was adjusted to match the asymptotic core collapse rate 
with the results of higher order fluid dynamical models (A — 0.186, Louis 1990), we chose 
here a value of A = 0.4977 (equivalent to C = 0.104), which yields fair agreement with 
Fokker-Planck results (Heggie & Stephenson 1988) as well as with our AT-body results 
in pre- and post-collapse. Throughout most of our calculations we chose Aa = 0.1, for 
reasons which will become clear in Sect. 4; a few runs were performed to show the effect 
of choosing Aa = 1- 

For our comparisons presented here we use only rather low particle numbers (A = 250, 500, 
1000, 2000); the post-collapse solution in such cases is a steadily reexpanding system with 
a regular core, related to the stable post-collapse isotropic self-similar models of Goodman 
(1987). A survey of larger particle numbers with the anisotropic gaseous model presented 
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here where the post-coUapse model is unstable to gravothermal oscillations is in progress 
(Spurzem & Louis 1993). Occasionally we will refer to another comparison of this model 
with an A'" = 10000 direct A'^-body calculation (Spurzem & Aarseth 1993). 

3 iV-BODY CALCULATIONS 

Over the last few years a new kind of computer became available for scientific use, 
consisting of a set of parallel transputers or vector processors. Their use opens a new 
chapter for stellar dynamical study. The power of these computers can be exploited in two 
different ways. One direction is simply to carry out calculations more quickly. This means 
to study larger systems. A second way is to obtain results having better statistics. It is 
well known that errors of the positions and velocities of the stars grow exponentially on 
a time scale shorter than the crossing time (Miller 1964, Goodman, Heggie & Hut 1993). 
Nevertheless A^-body modeling practitioners believe that statistical results obtained from 
N-body calculations are meaningful (cf. e.g. Aarseth & Lecar 1975). Therefore it is worth 
paying much more attention to improving the statistical quality of such models. This aim 
was achieved by running the same model several times in parallel, where the individual 
processes differ only with respect to the randomly generated initial positions and velocities 
of the stars. A full discussion of the method and the hardware used in the simulations of 
evolution of small A^-body systems (A^ = 250, 500, 1000, 2000) is described in great detail 
in the separate paper (Giersz & Heggie 1993a). Here we only outline them. 

The computational work has been carried out on two different parallel computers installed 
at the Edinburgh Parallel Computing Centre. One is the Meiko transputer array consisting 
of 400 processors, the other is the Meiko vector processor array containing 64 processors. 
Models consisting of A^ = 250, 500 and 2000 particles were computed on the transputer 
array using 56, 56 and 16 processors, respectively. The model for A^ = 1000 bodies was 
computed on the vector array using 40 processors, and on the DEC Alpha superclaster 
machine (called "fringe") for 20 cases. 

The first sets of models (A^ = 250, 500, 1000) were computed using a code developed 
mainly by Heggie (Heggie 1973). The A^ = 2000 model, not fully completed yet, and 
part of the A^ = 1000 model were computed using a version of Aarseth's standard code 
NB0DY5 (Aarseth 1985), as well as the N = 10000 case, which will be presented in more 
detail elsewhere (Spurzem & Aarseth 1993). The initial conditions for these models were 
drawn from Plummer's model (isolated system), with equal masses and no primordial 
binaries. In intervals of one A^-body time unit there was output produced for each of 
a parallel set of cases. This output mainly consisted of the following information: 1) - 
Lagrangian radii (1%, 2%, 5%, 10%, 20%, 50%, 75%) with respect to the centre of density, 
2) - the mean square radial and tangential velocities for each Lagrangian shell, 3) - global 
information: number of bound stars, energy of the system, number and energy of escapers, 
4) - information about binaries: radius, internal energy, number and energy of escapers. 
The output from all models (different processors) was collected in a single file and later 
analyzed (averaged over all models) to produce statistics for all important quantities at 
each time. 
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The improvement in the statistical quahty of the data achieved in our simulations allows 
us to perform very detailed comparisons with other methods, as in this case with the 
anisotropic gaseous model. 

4 RESULTS 

Three different groups of complete AT-body runs, consisting of AT = 250, 500, and 
1000 particles, have been compared with the anisotropic gaseous models; two not fully 
completed A^-body models will also be partly discussed (A^ = 2000, A^ = 10000). In what 
follows a phrase such as "1000-body model" refers to values averaged over many models, 
except for N = 10000 where there is only one case yet available. 

We compare Lagrangian radii containing certain fractions of the total mass as a function 
of time and the anisotropy of the velocity distributions averaged within these Lagrangian 
mass shells as a function of time for both the gaseous and A^-body models. All A^-body 
data had to be adjusted in their initial Lagrangian radii since the construction of an initial 
AT-body model, its scaling to the required value of the initial potential energy, and the 
use of the density centre all result in the fact that the initial model is slightly biased 
with respect to Plummer's model; the gaseous model was initially much closer to it. Fig. 
4 demonstrates the amount of shift necessary for the A^ = 1000 particle number as an 
example. 

Fig. 5 depicts the evolution of the 1% Lagrangian radius in the A^ = 1000 model in 
comparison to gaseous model results with a varying initial time tbo for the binary energy 
generation. The value of tbo controls the time at which the gaseous model curve deviates 
from a pure core collapse case without any energy generation. It does not alter very much 
the post-collapse behaviour, but we selected for the three particle numbers 250, 500, and 
1000 values of tbo — 50, 130, and 230, respectively, which appeared as best fits of the 
gaseous model curves to the A'"-body case. Note that the best fit is usually meant with 
respect to the 5%, 2% and 1% Lagrangian mass shells. We have chosen those radii because 
they are only slightly smaller than the core radius, within which most of the binary activity 
occurs. The evolution of the outer shells does not depend strongly on tbo- 

Physically tbo should be related to a time at which binary activity begins to play a role in 
the A^-body system; the ratios of the total binary energy released in the N-body system 
until the time tbo to the initial total energy of the system, however, were not the same 
for aU particle numbers (the fractions were 2.48%, 2.05%, 0.13% for N = 250, 500, 1000, 
respectively, see discussion of values below and in Sect. 5). 

Another parameter is the strength of the binary energy generation C^; Fig. 6a shows 
that the minimum of the curve for the innermost Lagrangian radius (1%) is clearly of a 
different shape for varying C^-values; for smaller Ci, the minimum is less shallow, because 
core collapse can proceed further to higher densities, until the energy generation suffices 
to halt and reverse it. There is also a variation in the total energy input for different 
Cb which can be recognized by the different post-collapse curves. Again for each particle 
number we could select those values of Cb = 55, 70, and 90, for N = 250, 500, 1000, 
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respectively, which appeared to fit best the shape of the minimum of the AT-body result 
as well as its post-collapse slope. Since there are always some fluctuations in the A^-body 
system we varied Ch only in larger steps (45, 55, 70, 90, 135) without attempting to 
determine the value more exactly. Moreover, one should note that due to a loss of cases 
in the parallel A'"-body calculations (for some individual models the total energy was not 
preserved due to very strong interactions between binaries and field stars and between 
binaries themselves (Giersz & Heggie 1993a)) the statistics of the corresponding averaged 
model worsen towards the end of the presented post-collapse evolution. For = 250, 500, 
1000 at the time t = 600 up to half of the initial number of cases were lost. 

According to Goodman (1987) (see his Eq. 11.14) Cb measures the product of the quantities 
X = 30c/^c5 where (f)c and Vc are the central potential and 3-D velocity dispersion, and /s, 
which describes the amount of energy supplied to the core by an individual binary. 

Fig. 6b depicts the evolution of x for difi'erent AT-body models. The time for all models 
N ^ 1000 was rescaled such that the distant two-body encounter relaxation timescale has 
the same value as for A'" = 1000; (this leads as one should expect to very good agreement 
of all pre-coUapse curves) . The curves were obtained by smoothing the original data using 
the standard procedure SMOOFT from Numerical Recipes (Press et al. 1986). 

Note that we find for A^ = 1000 a value of = 90, which is expected if the average total 
energy available from a binary resident in the core is supplied to the core (/a ~ 1). For 
lower particle numbers the changes of x with N during the post-collapse phase cannot 
be only accountable for the estimated values of Cb- This already suggests that for these 
models only a fraction, /a, of the energy released by 3-body binaries is deposited into the 
core, /a ^ 0.9 for A^ = 500 and /s fti 0.8 for N = 250, although we will give some more 
evidence on this point later. 

The best values for tbo and Cb were used to produce the data presented in Fig. 7 for the 
innermost Lagrangian radius (A^ = 250: Cb — 55, tbo — 50; A^ = 500: Cb — 70, tbo — 130; 
N = 1000: Cb = 90, tbo = 230), and in Figs. Sab to lOab for all three particle numbers 
showing the evolution of seven Lagrangian radii ranging from 1% to 75%. Looking at the 
innermost Lagrangian mass shells containing 1%, 2%, and 5% of the total mass we conclude 
that there is an excellent agreement between both models. Such a fitting procedure with 
the anisotropic gaseous model as described above would not be possible if the A^-body 
results were not improved in their statistics by computing multiple cases in parallel, which 
suppresses the intrinsic noise of the individual A'"-body system. 

Note, however, that for N = 250 the best fit to the post-collapse expansion made by 
C5 = 55 for the inner Lagrangian mass shells did not automatically produce the correct 
form of the minimum at core bounce in contrast to the other particle numbers A^ = 500 
and A^ = 1000; for A^ = 250 the minimum of the A'^-body model curves is always less 
shallow than for the gaseous model. 

For intermediate Lagrangian radii containing 10%, 20%, and 50% of the total mass at 
N = 1000 there is a tendency for the gaseous model to expand faster in post-collapse, 
at N = 500 there is fair agreement between both models here, whereas for N = 250 on 
the contrary the A'"-body models expand faster. The 75% radius is in all A'"-body models 
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further outside than in the gaseous model. The reasons for these remaining differences 
between our two models will be discussed also in Sect. 5. The larger expansion rate of the 
N-body system in the outer halo (75 %) as compared with the gaseous model is connected 
both with substantial amounts of escapers (Goodman 1984) and with degradation of the 
statistics due to loss of cases. There is a different treatment of the outer boundary in 
both models; the A'"-body model had initially all particles within a sphere of radius 10 
Th {rh' half-mass radius) i.e. any formed in the initial Plummer model at larger radii 
were removed, in the course of the evolution, however, a significant fraction of very loosely 
bound halo particles and escapers move outwards and populate areas within a sphere of a 
radius much larger than the initial radius; this process starts from the very beginning with 
particles just moving outwards on radial orbits. The phase space is not completely covered 
in that region, e.g. the angular momentum of the particles there (which are approximately 
5% of all particles for N = 10000 in the late collapse phase) is always much less than 
the maximum angular momentum of a circular orbit. On the contrary the gaseous model 
extended much further outwards already in the initial model in order to avoid dynamical 
perturbations originating from the boundary, but the system was assumed to be isotropic 
there initially, i.e. the accessible phase space was fully covered. Any expansion of the 
gaseous model's 75% radius is a dynamical expansion driven by pressure against the outer 
shells, which is physically different to what happens in the A^-body system, where some 
fraction of stars is just moving outwards filling empty space. Therefore the 75 % Lagrangian 
radius is always further outside in the A^-body models. This effect is even more pronounced 
for the 90% radius of the A^ = 2000 body calculation as an example (Fig. 11), where one 
can recognize that the A^-body model started with a smaller radius than the gaseous model, 
but expands from the beginning faster, since particles can move freely outwards in contrast 
to the gaseous model. The 90 % radius has therefore not been shifted, since its difference 
with the gaseous model was much larger than usual. 

We now start to discuss the evolution of the anisotropy A = 2 — 2cr^/cr^; it is presented 
as a mass- weighted average taken over the spherical shells whose inner and outer radii are 
given by the Lagrangian radii of the previous figures. It is now the only remaining free 
parameter Aa which limits the anisotropy by coUisional isotropization. Let us denote with 
As, Aq, A'j the anisotropy between Lagrangian radii of 10% to 20%, 20% and 50%, and 
50% to 100%, respectively. Their time evolution is shown in comparison between N-hody 
results and gaseous models for N = 250, 500, 1000, 2000, and 10000 in Figs. (12) to (16), 
respectively. Most gaseous models were run with Aa = 0.1 (solid lines in Figs.), whereas 
in Figs. (12)-(14) alternative results for Aa = 1 (dashed lines) are also presented. The 
anisotropy generation for the Aa = 1 case is clearly much too high. Aa = 0.1, however, 
yields excellent agreement of ^5, Aq and A^ for N = 10000 (Fig. 16), as well as for all 
other particle numbers until a certain time tao, which turns out to be approximately equal 
to tbo for A^ = 1000, and a little smaller than tbo for A^ = 500 and A^ = 250. Since for 
the lower particle numbers A^ < 10000 we know that there has been already more binary 
activity at tbo it is consistent to postulate a connection between the binary activity in the 
core and the larger anisotropy (more radial energy) outside. 

There is further evidence that the anisotropy in the outer shells is determined (at least 
partially) by binary activity. Firstly, there is a strong increase in the anisotropy after tbo, 
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when the binary energy generation starts, in all Aa = 1 gaseous models (Figs. 12-14). For 
the Aa = 0.1 case it is suppressed by stronger coUisional isotropization. Secondly, the total 
energy and number of the bound hard binaries level off during the post-collapse evolution 
nearly at the same time as anisotropics do (Giersz & Heggie 1993b). 

To further clarify the point we performed an artificial experiment by substituting the 
density exponent of 3 in Eq. (9) with 3 — /9 for the tangential and 3 + (3 for the radial 
energy generation. We compared runs for /? = 1 and /? = 0, the latter being the standard 
case. The total energy generation was normalized such that it does not differ in both cases 
and such that it is still isotropic at the centre. For positive f3 there is more radial binary 
energy generation than tangential; the physical picture is that encounters with binaries 
very likely produce a particle leaving the core with low angular momentum i.e. higher 
radial than tangential energy. It can be seen in Figs. (17) and (18) that this indeed is a 
good model for what happens with the anisotropy: the stronger radial energy generation 
does not alter and Aq much, but there is a clear trend towards the A'"-body results for 
Ar. 

As a last check we suppressed close encounter and binary activity in an A^-body simulation 
using S. Aarscth's NBODYl with a smoothing parameter e = po and e = 2po, where po is 
the impact parameter for which the deflection in the relative orbit of interacting particles 
is equal to 7r/2. In A-body units po = 2/(AF^), where for V'^ we use a value obtained 
from the virial theorem. Assuming that the maximum impact parameter is roughly equal 
to Th we can estimate for N = 250 that for e = 0.016 and e = 0.032 about 9% and 21% of 
all encounters is suppressed, respectively. 

The striking result shown in Fig. (19) is the much better agreement between A-body and 
gaseous models for A^ in the case e — 0.016; the trend is a little too strong for e = 0.032. 
Taking this last result into account we conclude that Aa = 0.1 seems to be the physically 
realistic value for all cases, provided close encounters and binary activity play little role; 
this is naturally the case in large A systems (for which the gaseous models are qualified). 
It is consistent that in all results for Aa = 0.1 and t < t^o fhe agreement is excellent; any 
difference starts for low A and at times after which binary activity has started. 

5 CONCLUSIONS AND DISCUSSION 

We have performed a set of direct A-body models of idealized single mass star clus- 
ters with improved statistics obtained by calculating independent models in parallel on a 
parallel computer for particle numbers of A = 250, 500, 1000, and 2000. One model of 
A = 10000 (Spurzem & Aarseth 1993) complements some of our data. The results of the 
parallel A-body project are published in more detail in Giersz & Heggie (1993ab); here 
we related them to time dependent evolutionary calculations with an anisotropic gaseous 
model, which is closely related to the models of Bettwieser & Spurzem (1986) and Louis 

6 Spurzem (1991, LS), but complemented by an energy source appropriate to describe a 
(usually isotropic) energy input by formation and hardening of three-body binaries analo- 
gous to that used in the isotropic gaseous models of Bettwieser and Sugimoto (1984) and 
Heggie and Ramamani (1989). 
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Excellent agreement of the evolution of the innermost Lagrangian mass shells (less than 
or equal to the core radius) in both pre- and post-collapse could be found by adjusting 
the strength and start time of the binary activity in the gaseous model. Remarkably the 
numerical value of the binary strength Cf, = 90 found in that manner agrees very well 
for N > 1000 with theoretical expectations for systems whose evolution is dominated by 
small angle gravitational encounters (Goodman 1987). For = 500 and = 250 such 
agreement could be achieved only by choosing somewhat smaller values of = 70 and 
Cb = 55, respectively. 

The generation of anisotropy during pre-coUapse and its levelling off in post-collapse within 
the half-mass radius is in fair agreement between A^-body and gaseous models, provided 
one takes as timescale for collisional decay of anisotropy XaTa, where Aa = 0.1, and Ta is 
the local anisotropy decay time for an anisotropic Larson-type distribution function (see 
Eq. (5), (6), and (A15), and Larson 1970). For the outer portions of the system the 
agreement is satisfactory only if N > 2000 or time t < tbo, the latter being the time at 
which binary activity becomes important. 

It is important to note that in pre-coUapse and without binary activity we have scaled 
all A^-body and gaseous models such that their collapse phases match exactly; the scaling 
factor was computed theoretically and the result is that all models for all particle numbers 
agreed excellently for a unique value of the remaining parameters (A = 0.4977 for heat 
conductivity, A a = 0.1 for anisotropy decay, and 7 = 0.11 in the Coulomb logarithm, 
cf. Eqs. 2-8). The agreement shows up for the anisotropy everywhere in the system 
and for the evolution of the Lagrangian mass shells up to the core radius. It is a puzzle, 
that between the core and half-mass radius in pre- and post-collapse the evolution of the 
Lagrangian mass shells in the gaseous model does not agree very well with the AT-body 
models. There is a tendency for the gaseous model to expand faster (this is also the case 
for the isotropic Fokker-Planck model, Giersz & Heggie 1993b) at these radii as can be seen 
for the A^ = 1000 case in Fig. (8ab); however the same result occurs in the runs available 
for the largest particle numbers A^ = 2000 and 10000 (no figures given in this paper). 
Especially for N = 10000 distant two-body encounters are the main force which drive 
the evolution. So the gaseous model based on the thermal conductivity approximation 
should in principle describe the evolution well. However, this is not the case, particularly 
for the 50% Lagrangian radius. An explanation can be connected with the fact that the 
gaseous model cannot cope with non-local scattering. To mimic the energy transport in 
the radial and tangential directions we have to take into account all interactions along a 
particle orbit. This could be done for example by an orbit- averaged 2-D Fokker-Planck 
equation or by Monte-Carlo realization of this equation. We would like to stress that the 
best agreement with A^-body data for the half mass radius is achieved by Stodolkiewicz's 
Monte-Carlo code (1982). Unfortunately any quantitative proof of this statement is not 
possible because the only obtainable data are diagrammatic. 

For smaller systems (A^ = 250 and 500) that effect is present as well but we cannot see 
it because it is gradually overcome by the effect of reaction products of close two- and 
three-body interactions, which deposit their energy by large angle encounters in the outer 
regions of the system. 
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One may argue that a variation of 7, the factor in the Coulomb logarithm, or A, the heat 
conductivity coefficient, should be checked to see whether it is sufficient for full agreement 
of the evolution of all Lagrangian mass shells. There is some theoretical evidence about 
a possible dependence of 7 on position and the total number of stars (see Spitzer 1987 
for detailed discussion). Spitzer (1987) suggested that 7 could be equal to 2Nc/N. From 
Heggie & Stephenson (1988) one may conclude that A is slightly different between pre- 
and post-collapse. However, our results indicate that our adopted constants give very good 
agreement during the collapse phase. So the local effects of distant two-body encounters are 
correctly modeled in this phase. Discrepancies start to build up when the number of stars 
in the core is very small and binary activity starts to play a role. This can favour Spitzer's 
definition of 7, but on the other hand non-local effects connected with the small number 
of stars start to be increasingly important. For example very close two- and three-body 
scatterings can produce high velocity stars which would disturb the velocity distribution 
such that our second order moment equations are not sufficient to model the system. Our 
anisotropic gaseous model improves previous gaseous models with respect to the inclusion 
of (second order) anisotropy, but it may be too simplified in the presence of high-energy 
scattering events, which happen in the core predominantly for small particle numbers 
< 1000 in post-collapse. It would be interesting to check how an anisotropic gaseous 
model such as that of Louis (1990) or numerical solutions of the orbit-averaged 2-D Fokker- 
Planck equation could cope with this problem. So we suggest that post-collapse evolution 
in principle can be understood without introducing new values of 7 and A, although we 
cannot totally rule out that possibility (particularly for high N systems) on the basis of 
our results. 

There is more evidence that binary and close encounter activity causes differences between 
gaseous and A^-body models for low N; this is the much too strong anisotropy generation 
outside of the half-mass radius in these cases (Figs. 12-14); first the problems occur as 
the binary activity and the anisotropy levels off in the A'"-body systems at the same time 
as the binary activity begins (Giersz & Heggie 1993b). Second, a test calculation with a 
local, but anisotropic (stronger radial) energy generation due to binaries exhibits a trend 
towards the A^-body results. Third, we may interpret our empirically found strengths of 
binary energy generation in terms of the parameter /s of Goodman (1987) describing the 
fraction of binary energy liberated in the core. Our data indicate that /s < 1 for A?" < 500. 
Since all binaries are located in the core this means a rapid non-local energy transfer to 
regions outside of the core radius, carried by particles moving outwards on radial orbits 
and creating the anisotropy in the halo. Fourth, a test N-body calculation with suppressed 
close encounters by choosing non zero smoothing parameter shows much better agreement 
with the gaseous model results. 

Our findings are consistent with the results of earlier work of Spitzer & Mathieu (1980) 
and Goodman (1984) stating that the reaction products of superelastic binary-single star 
scatterings will be transported outwards on elongated radial orbits, sometimes even escape. 
Angular momentum scattering along these orbits will distribute the energy of the reaction 
product in the corresponding radial zones. This is the mechanism by which the above 
mentioned non-local energy transport is provided and we identify the reaction products as 
carriers of the observed high anisotropy in the A'"-body models. It is not surprising that 
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for systems with larger N and higher x (deeper central potential, see Fig. 6b) this energy- 
distribution is confined more to the homogeneous core itself and thus more consistent with 
the gaseous model. 

The clue to an understanding of the remaining further differences between the two models 
(sharp minimum of the Lagrangian radii, compare Figs. 7 and 10a, and higher binary 
energy generation for small N) lies in the nature of the binary and close encounter activity, 
whose character is very different for the low N systems and for high N systems. For 
the former there is a higher probability, with respect to the probability of small angle 
encounters, that core particles are subject to a close two-body encounter or a three-body 
binary will form. Therefore for = 250 binaries are created earlier and they can affect 
the core collapse and eventually stop it earlier (in terms of x = 3(f)o/vQ, see Fig. 6b) than 
in higher N systems. Because of smaller x the rate of collapse is higher (Cohn 1980), 
and binaries have to generate more energy to influence core collapse and finally reverse it. 
Thereby we finally have found a natural explanation why the energy generated by binaries 
before the time t^o is higher. However, due to smaller x (less deep central potential) a 
bigger fraction of the energy can be carried by particles non-locally outwards and it is 
not felt as early by the core as for higher particle numbers. This explains also why the 
intermediate shells for N = 250 expand faster than in the gaseous model, because that is 
where much of the energy and mass is non-locally transported to. 

A supplementary explanation for the sharper minimum in the time evolution of the inner 
Lagrangian radii in the N-body models (see Fig. 10a) is connected with the growing 
importance of statistical fluctuations with decreasing number of stars. For N = 250 the 
number of stars in the core at the time of core bounce is extremely small (about 14 
particles). Therefore stochastic processes connected with formation of binaries and their 
subsequent burning cannot be properly approximated by continuous formulae (Goodman 
1984, 1987). Results obtained by Giersz & Heggie (1993b) for stochastic binary formation 
and burning with an isotropic gaseous model suggest that the minimum is in this case 
much sharper than in the standard case, whereas the post-collapse expansion is nearly the 
same. On the basis of that result we conclude that statistical fluctuations are important 
and at least partially account for the sharper minimum. 

It can be concluded that the core evolution for large particle numbers is excellently de- 
scribed by anisotropic gaseous models in pre- as well as in post-collapse with the stan- 
dard phenomenological implementation of the binary energy fed into the core. Remaining 
problems are faced by the gaseous model in two respects: i) small particle numbers and 
three-body encounters: they tend to spread energy non-locally accompanied by a large gen- 
eration of anisotropy, which cannot properly be taken into account by the present simple 
isotropic energy input formula for binaries in the gaseous model; ii) large particle numbers 
and two-body scatterings: in the regions outside the core of the system an appreciable dif- 
ference occurs in the rate of evolution of the Lagrangian radii - the gaseous model evolves 
too fast. Whether this is a result of approximations inherent in the gaseous model (e.g. 
local approximation for all coUisional effects, or neglect of higher order moments of the 
velocity distribution) remains a question for future work. One could for example take into 
account non-local scattering in the context of a numerical solution of the orbit-averaged 
2-D Fokker-Planck equation or a more detailed form of the velocity distribution by a higher 
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order moment model. 
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APPENDIX 



Derivation of the collisional decay of anisotropy 

According to Larson (1970) we approximate the real local (r-dependent) velocity 
distribution function / as a function of the modulus of a star's random velocity vector 
V and the cosine of its angle with the radial direction /i, 

oo 

i=o 

where Pj^jj) denotes the Legendre polynomials of order j. We are going to apply the 
Fokker-Planck equation in v, ji coordinates, which was computed by Rosenbluth, McDonald 
& Judd (1957, RMB); therefore we need to know the Rosenbluth potentials 

/i(v) = 2m j /(v')|v - Y'\-^d\' (A2) 

^(v) j /(v')|v - v'|dV (A3) 
for the distribution function of Eq. (Al); they are 

oo 

h{v, I^)-S^m.f2^ [ ^Mv'WV (A4) 
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g{v,fj,) =^^rn-^——j- 

1=0 

oo 
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oo > 



(A5) 







with v< = min(v,v'), v> = max(v,v'). 

Since we are interested mainly in distribution functions with anisotropy in the second order 
moments the series in Eq. (Al) shall be truncated now at that order and we use in accord 
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with Larson's (1970) definitions 



^0 = giv) 
^1 = 



with 



C2 = 



^2 ^2 



3(72 • 



{A7) 



is the isotropic Maxwell-Boltzmann distribution. These definitions just ensure that the 
quantities cr^, cr^, and 3a^{vr- — u) are recovered by determining the second and third order 
moments of /; fourth order moments have been assumed to take the same value as in the 
case of a Maxwell-Boltzmann distribution function; in Larson's notation this means that 

we have set ^ = (see his Eq. 2). 

Now the Rosenbluth potentials can be further evaluated in terms of the functions 
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As a final result we get 

h{v, ji) = 87rm 

g{v, n) = 47rm< 
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With these potentials the Fokker-Planck equation of RMB is determined; let the right 
hand side of their Eq. (31) be denoted with FP; we then evaluate 
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It turns finally out that 

n n Pr-PtfQ Q Pr-Pt\ /...x 

with the average pressure p = {pj. + 2pt)/'i = pa^. Since the second term in the above 
equation is small for reasonable values of the anisotropy we finally get 

Dpr-Dpt 9 ^ ^ 



22 



FIGURE CAPTIONS 



Fig.l Evolution of Lagrangian radii containing the stated percentage of the total mass for 
an isotropic (IGM) and anisotropic (AGM) gaseous model for a N = 1000 star cluster 

Fig. 2a Radial profile of the logarithmic density power law index a for a number of gaseous 
models with different evolutionary time. The dash-dotted line represents the self-similar 
anisotropic pre-coUapse model of LS 

fig. 2b Radial profile of the anisotropy A for a number of gaseous models with different 
evolutionary time. The dash-dotted line represents the self-similar anisotropic pre-coUapse 
model of LS 

Fig.Sa Comparison between the density power law index a as a function of AaA for 
anisotropic gaseous models (AGM) and self-similar anisotropic pre-collapse models of LS 
(SIM). The values of a for the isotropic gaseous model of Lynden-Bell and Eggleton (LBE) 
and the isotropic Fokker-Planck model of Cohn (FP) are marked 

Fig. 3b Comparison between the anisotropy A as a function of A^A for anisotropic gaseous 
models (AGM) and self-similar anisotropic pre-collapse models of LS (SIM) 

Fig.4 Amount of shift necessary to adjust the averaged N = 1000 A^-body model to the 
anisotropic gaseous model for 1% Lagrangian radius 

Fig.5 Evolution of the 1% Lagrangian radius in the averaged N = 1000 A'"-body model in 
comparison to the anisotropic gaseous model for different initial times tbo for switching on 
the binary energy generation. The cc indicates a pure core collapse gaseous model 

Fig. 6a Evolution of the 1% Lagrangian radius in the averaged N — 1000 A'"-body model 
in comparison to the anisotropic gaseous model for different strength of the binary energy 
generation parameter Cb- The cc indicates a pure core collapse gaseous model 

Fig. 6b Time variation of x (Def. see text) for three particle numbers indicated at the 
curves, in the averaged A'"-body models; all times were normalized such that the core- 
collapse time is the same for all N and the strongly fiuctuating curves were smoothed. 

Fig. 7 Evolution of the 1% Lagrangian radius in the averaged N = 1000, N = 500, N = 250 
A^-body models in comparison to the anisotropic gaseous models (C5 = 90, ti,o = 230), 
{Cb = 70, tbo = 130), {Cb = 55, tbo — 50), respectively 

Fig.Sa Evolution of the 1%, 2%, 5%, 10% Lagrangian radii in the averaged N = 1000 
A^-body model in comparison to the anisotropic gaseous model 

Fig.Sb Evolution of the 20%, 50%, 75% Lagrangian radii in the averaged N — 1000 
AT-body model in comparison to the anisotropic gaseous model 

Fig.9a The same as in Fig. 8a but for N = 500 

Fig.9b The same as in Fig. 8b but for N = 500 

Fig. 10a The same as in Fig. 8a but for N = 250 

Fig.lOb The same as in Fig. 8b but for N = 250. The 75% Lagrangian radius is not 
available for this model. 

Fig. 11 Evolution of the 90% Lagrangian radius in the averaged N = 2000 AT-body model 
in comparison to the anisotropic gaseous model 
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Fig.l2 Evolution of the anisotropy between Lagrangian radii of 10% to 20% (A5), 20% to 
50% (Ag), and 50% to 100% (A7) in the averaged = 1000 A^-body model in comparison 
to the anisotropic gaseous models for different parameters of coUisional anisotropy decay, 
Xa 

Fig. 13 The same as in Fig. 12 but for N = 500 
Fig. 14 The same as in Fig. 12 but for = 250 

Fig. 15 Evolution of the anisotropy between Lagrangian radii of 40% to 50%, 50% to 75% 
and 75% to 90% in the averaged N = 2000 A/"-body model in comparison to the anisotropic 
gaseous model without binaries; note the different definition of Lagrangian shells here. 

Fig. 16 The same as in Fig. 12 but for N = 10000 and Xa = 0.1 and no binaries 

Fig. 17 Evolution of the anisotropy between Lagrangian radii of 10% to 20%, 20% to 
50% and 50% to 100% in the averaged N = 1000 A'"-body model in comparison to the 
anisotropic gaseous models for different /3 (see text) 

Fig.18 The same as in Fig. 17 but for = 500 

Fig. 19 Evolution of the anisotropy between Lagrangian radii of 50% to 100% in the 
averaged A^ = 1000 A'^-body models with different smoothing parameters e in comparison 
to the anisotropic gaseous model 
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